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Abstract. We propose and analyze a goal-oriented a posteriori error estimator for the atomistic- 
continuum modeling error in the quasicontinuum method. Based on this error estimator, we develop 
an algorithm which adaptively determines the atomistic and continuum regions to compute a quan- 
tity of interest to within a given tolerance. We apply the algorithm to the computation of the 
structure of a crystallographic defect described by a Frenkel-Kontorova model and present the re- 
sults of numerical experiments. The numerical results show that our method gives an efficient 
estimate of the error and a nearly optimal atomistic-continuum modeling strategy. 



1. Introduction 



The quasicontinuum (QC) method [22, l23|, |2j] has been successfuhy used to efficiently cou- 
ple atomistic and continuum models for crystalline solids and offers the possibility of computing 
mesoscale or macroscale properties by a nearly minimal number of degrees of freedom. Accurate 
modeling requires that an atomistic model be used in regions with highly non-uniform deformations 
such as around dislocations, whereas a continuum model can be used in regions with nearly uniform 
deformations to reduce the number of degrees of freedom. 

It is usually not known a priori which regions of some specimen undergo uniform deformations 
and which do not, so a posteriori error estimation is important for the design of efficient numerical 
approximations by the quasicontinuum method. Since the purpose of a computation is often to 
obtain the value of a (usually local) quantity of interest to a desired error tolerance rather than to 
obtain a solution to a desired error tolerance for a global norm, there has been great interest in 
the development of goal-oriented error estimators for many problems. They are based on duality 
techniques and have been developed and used to adaptively refine finite element approximations of 
continuum problems [l|, 0] and to study and control modeling error [l^ . 

In this paper, we extend this approach to develop an a posteriori error estimator for the qua- 
sicontinuum method which quantifies the atomistic-continuum modeling error for a goal function 
and allows for an adaptive decision about which regions can be accurately modeled as a continuum 
and which regions need to be modeled atomistically. Methods to determine the optimal mesh size 
within the continuum region will be studied in a forthcoming paper. 

Crystallographic defects provide a challenge to validate atomistic-continuum error estimators 
and adaptivity. No such error estimators and adaptive methods currently exist for fully three- 
dimensional crystals. As a step in this direction, we develop a rigorous theory for a simple one- 



dimensional atomistic model for a defect that is a modification of the Frenkel-Kontorova model [151 ]. 
We add next-nearest-neighbor harmonic interactions between the atoms to the nearest-neighbor 
harmonic interactions between the atoms in the classical Frenkel-Kontorova model. 

A priori analyses for various quasicontinuum approximations have been given in 0, [lol . 



121 . Il3l . [ij, |2l|]. An a posteriori analysis for a slightly different one-dimensional quasicontinuum 
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approximation is given in [20[ • The development and application of a goal-oriented error estimator 
for mesh coarsening in a two-dimensional quasicontinuum method is reported in 17|, ll8|] . 

Let us mention that the continuum model used in the QC method, which coincides with the 
model obtained by the classical thermodynamic limit, is by far not the only reasonable continuum 
model to use. A method to derive continuum models which approximate atomistic models up to 
an arbitrarily high order has been proposed in Q]- 

The paper is organized as follows. In Section [2l we give a general formulation of the one- 
dimensional quasicontinuum approximation [1^ that includes not only two-body and three-body 
potentials, but also many body potentials such as the embedded atom potential [u, 7]. In Section[3l 
we describe our extension of the Frenkel-Kontorova model and its quasicontinuum approximation. 
In Section HI we introduce the primal and dual problems for our model and formulate our approach 
to goal-oriented error estimation. 

Next, in Section [5] we extend the approach in |16j] to develop an error estimator for atomistic- 
continuum modeling. This first error estimator does not allow a decomposition among the atoms 
that can be used for atomistic-continuum adaptivity, so we propose and analyze a less accurate 
second error estimator that does allow such a decomposition. 

Finally, in Section [6] we propose an adaptive atomistic-continuum modeling algorithm and show 
that it gives an efficient estimate of the modeling error and a nearly optimal atomistic-continuum 
modeling strategy for the computation of defect structure. 



2. Quasicontinuum Approximation 

The departure point for the QC approximation is the potential energy of the atomistic system. 
The potential energy that is utilized fully models the properties of the system. The local minima 
of the potential energy model the metastable states of the system, and the potential energy can be 
used in Newton's equations of motion to model the dynamical behavior. 

The QC method approximates the potential energy of the atomistic system in two steps. First, 
we develop a continuum potential energy that will be used in the adaptively determined continuum 
region, and we then show how to reduce the degrees of freedom in the continuum region. 

2.1. The Atomistic System. We assume that the atomistic system has 2M atoms with deforma- 
tion given by = {y°Lj^,ij^i, ■ ■ ■ y%i) € M^*^. Without loss of generality, we assume that the atoms 
are ordered so that their positions satisfy < y'^j^i- Furthermore, we assume that the atomistic 
total potential energy, f "(y"), can be written as a sum over potential energies associated with each 
atom, ^"(y"), so that 

M 

^"(y")= E ^"(y")- (2.1) 

i=-M+l 

This decomposition can be found for most empirical potentials, including embedded atom potential 
energies [1, 0]- For example, if the atomistic total potential energy £'"(y") is given by 

where ip{r) is an empirical two-body potential energy, then we can obtain the decomposition ()2.ip 
by taking 

w) = ^Ev'(y"-2/n- (2.3) 
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We note that 5f(y'^) can also contain contributions from external forces, such as for the Frenkel- 
Kontorova model described in Section [3l and can thus depend on i. 

2.2. The Atomistic-Continuum Energy. For any deformation y" G M^^, we let L'-^+iy'^ g 
denote the linear extrapolation of the atomistic positions yf and y^_^i given by 

(r.^+V")fc = ik- i)yt+i + (i + 1 - k)yf for = -oo, . . . , oo. (2.4) 

The continuum potential energy £'?(y") of atom i is obtained from the average of the atomistic 
potential energy £^ evaluated at the extrapolations L*~^'*y'^ and L^'^+^y"^ by 

Sfiy-) := ifr(L-i V) + i^r(L*'*+V"), (2.5) 

where we note that the domain of £f has been expanded to the infinite periodic atomistic systems 
in the range of L*~^'* and We assume that f ? is finite for infinite periodic atomistic systems, 

which is true for (j2.3p when the two-body potential ipir) decays fast enough so that Ylk^i '>P{kr) is 
finite for r 7^ 0. At the endpoints of the chain, the extrapolation can be done only to one side, so 
we neglect the undefined part and define 

S^i+iiy") ■■= \£lM+i{L~''^''-''^h^) and £h{y-) := \£l,{L"'~^^'^y-). (2.6) 

We then decide for each atom i whether to model its energy atomistically by £f{y°') or as a 
continuum by £^{y°'). We thus obtain for the whole chain the atomistic-continuum energy 

M M 

i=~M+l i=-M+l 

M M M-l ^ ' ' 

i=-A/+l i=-M+2 i=-M+l 

where 

I 1 if atom i is modeled atomistically, , ^, ^. , , 

6i = { and Sf = I - 6i. (2.8) 

I if atom i is modeled as continuum. 

This approximation allows for a slightly faster evaluation of the energy and its derivatives, especially 
if f ? is long-ranged. However, it reveals its full strength only after the quasicontinuum coarsening to 
be described next. We note that sometimes atomistic degrees of freedom and energies are referred 
to as nonlocal and continuum degrees of freedom and energies are referred to as local [2^ . 

2.3. Repatoms: Reduction of Degrees of Freedom. The quasicontinuum method allows a 
reduction of the number of degrees of freedom in the continuum region. To this end, we choose 
so-called representative atoms, or more briefly called repatoms. The repatoms are a subset of the 
original atoms. The quasicontinuum approximation of the energy is defined completely in terms of 
the repatoms. 

We choose the repatoms by defining indices ij for j = —N + 1, . . . , N where 

-M + 1 = £-N+i <■■■ <ij < ij+i < ■ ■ ■ < In = M. 

The atoms at for i = £-]\f+i,i-N+2, ■ ■ ■ An are repatoms, and all of the remaining atoms are 
non-repatoms. We have that 

v,=l,+^-l, (2.9) 
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gives the number of atomistic intervals between the repatoms ij and ij+i- We require that the chain 
not be coarsened in the atomistic regions, which precisely means that 5^ = 6p ,,=... = 6p =1 
whenever Uj > 1. 

Finally, the interactions of the atomistic energy only partially reach into the continuum part if 
the atomistic potential has a finite cutoff radius. To allow for an exact calculation of this energy 
without atomistic interpolation, we require that these regions are not coarsened as well. As we 
will see in the next subsection, the atomistic next-nearest-neighbor interactions from the Frenkel- 
Kontorova model studied in this paper reach two atoms into the continuum part. Hence, we require 
that i'j-2 = ^j-i = = T^j+i = 1 whenever (5°, = 1. Other potential energies in general require 
similar conditions that depend on their cut-off radius. 

We denote the position of the j-th repatom by = and the vector of all repatoms by 

2.4. The Quasicontinuum Energy. Now we define the quasicontinuum energy. To this end, the 
missing non-repatoms are implicitly reconstructed. We will see later that this helps to set up the 
QC model, but needs not be done for the actual computation. 

The reconstruction is done by a linear interpolation between the nearest repatom to the right 
and to the left. That is, the vector of all atomistic positions is computed from the vector y"^"^ of 
repatom positions by 

J : M^^ ^ M2*^ (/y^^).,+^ := + -y]\,, m = 0,...,u,. (2.10) 

We note that 

yf = {Iynir (2.11) 

The underlying idea is that in regions where the lattice spacing of the atoms is nearly constant, 
this interpolation is very close to the actual atomistic positions and therefore leads to a good 
approximation of the total energy. Only a few repatoms are needed in these regions. This exactly 
corresponds to mesh coarsening in classical finite element approximations of continuum models. On 
the other hand, in regions where the lattice spacing is non-uniform, such as around a dislocation, 
all atoms must be chosen to be repatoms to obtain sufficient accuracy. This guarantees that the 
full resolution of the atomistic model in the critical regions is retained and corresponds to a high 
refinement in classical finite element continuum models. 

We define the QC approximation of the total energy to be 

£<l^(y1^)~e<^^{Iy1^), (2.12) 

Now (]2.12p has to be reformulated such that it can be computed efficiently, without the overhead 
of the interpolation. Most atomistic potentials are invariant to translations, a property that al- 
lows us to simplify (j2.12p considerably. For any translationally invariant energy 8f, we have that 
5f (L^'^+V) = (t>i{y'}+i - Vi) and (L^-^'V) = Mvf - vt-i) for some function 0^. If these 
functions (pi coincide, that is, (pi = for all i and we can write 

ga^^L^-l,^a^ = ^ {^y- _ y-_^) and (L^'^+V") = 4> fei - yt) (2.13) 

for some function : M — > M. Here (p plays the role of a continuum energy density and is given for 
the two-body potential (|2.2p by 

oo 
k=l 
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Equations ^J^, ([2J2]) . and ([2J3]) lead to 

M M M-l 

i=-M+l i=-M+2 i=-M+l 

M M 

= E miiyn + h E miiy'li - iiy'%-i) (2.14) 

i=-M+l i=-M+2 
M-l 

+ i E miiy'^'Ui-iiy'^h)- 

i=-M+l 

Because ly'^'^ is the linear interpolation between two repatoms and y'j\i, we have 

qc qc 

{ly^r+i - ily'^h = ^^±i^, i = . . - 1. (2.15) 



Hence, 



M N-l / qc _ gc\ 

£'\yn= E ^t^t{iy"')+ E ^^'^ V ^^-^^^ 

i=-A/+l j=-Af+l V / 



with weight factors 



if both yf^ and y^li are atomistic, 



'^i ~ + '^^j+i) ~ 1 I exactly one of y^ and y^+i is continuum, (2.17) 

i^'j if both y^ and are continuum. 

The first sum corresponds to the atomistic region which will be a small region and is thus com- 
putationally inexpensive. The second sum only involves at most 2N terms which is a considerable 
reduction when N <C M. 

Note that the second term in formula ()2.16p coincides with an integral over the energy density (j) 
as it occurs in finite element discretizations of classical continuum mechanical models. Hence the 
apparently unmotivated definitions (12. 7p and (j2.5l) of the continuum energy here result in what is 
commonly understood as a continuum energy. The linear interpolation operator / resembles the 
Cauchy-Born hypothesis. 

3. Frenkel-Kontorova Model 

Dislocations are lines in crystals which represent a defect in the lattice structure [ill, see Fig- 
ure [H Typically, there is a core of small radius surrounding the dislocation line where the lattice 
structure is highly deformed, but the lattice structure is nearly uniform outside the core. A sim- 
ple one-dimensional model for a defect such as a dislocation is given by the Frenkel-Kontorova 



model 15|. Here, the elastic energy is modeled by harmonic interactions between the atoms in the 
one-dimensional chain and the misfit energy of the slip plane is modeled by a periodic potential. A 
more accurate model of the same form is given by the Peierls-Nabarro model 

3.1. Atomistic Frenkel-Kontorova Model. We study a single defect in the middle of the chain 
of 2M atoms. To achieve a symmetric description in terms of bonds, we number the atoms from 
—M + 1 to M. The defect is situated between the atoms numbered and 1 (Figure [2]). 

Recall that the atomistic positions are denoted by = (y" . . . , y^) G M^*-'^. The total 
potential energy for this atomistic system is then a function E"" : M?'^ ^ M of the atomistic 
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Figure 1: Cross-section through a dislocation in a three-dimensional cubic lattice. The displayed plane 
repeats periodically in the three-dimensional crystal. Vertical bonds are shown by lines to emphasize the 
topological defect. 



o o o o 



O O O O 



-Ma„ 



■■ atomistic 



(-K-1)a„-Ka„ 
O = continuum 



Ka^ (K+1)aj, Ma^ 



Figure 2: Numbering of the atoms. The dislocation is situated in the middle of the chain between atoms 
and Di- 



positions. For the Frenkel-Kontorova model, the energy, f " = S""'^ + S""'"^, consists of two parts, 
namely the part which models the elastic energy of the defect, E""'*^ , and the part which models the 
misfit energy on the slip plane, 8°"'^. 

The elastic energy is modeled by Hookean (harmonic) springs between nearest-neighbors (NN) 
and next-nearest neighbors (NNN), and the total elastic energy is given by 

M-l M-1 

^a,e(ya)^ ^ ik,(y^^^-yf-aof+ J] ^,h{yt+, - y?^, - 2aof , (3.1) 

j=-M+l i=-M+2 

where the moduli ki > and k2 > describe the strength of the elastic interactions, and where 
qq £ denotes the equilibrium distance. 

We note that the asymptotic expansion to second order of any nonlinear NN/NNN potential 
energy 

Af-l M-l 



i=-M-|-l 



-M+2 
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Figure 3: Frenkel-Kontorova model. The wells depict the misfit energy 



about a" = [(-M + l)ao, (-M + 2)ao, • • • , (M - l)ao, Mao]^ € M^M has the form 



i=-M+l 

- iP'{2ao){yh - Vli-i - ao) - V''(2ao)(2/^A/+2 " V-M+i " oo) (3.3) 

Af-l A/-1 

+ l<(ao) 5] (y^^^_y«_ao)2 + iV'"(2ao) J] (yf+i - y^i " 2ao)2. 

We thus see that the elastic energy (j3.ip with ki = ip"{ao) and k2 = 'ip"{2ao) approximates the 
energy (j3.2p to second order if t{j'{ao) + 2-0' (2ao) = and if we ignore the boundary terms in the 
second line of (jS.Sp . 

The misfit energy of the shp plane is modeled by a periodic potential (Figure [3]). We model this 
misfit energy by 



J2 1^0 ( yf - ao 



^ + 1 

oo 2_ 



(3.4) 



where [xj denotes the largest integer smaller than or equal to x, and where the constant /cq > 
determines the strength of the misfit energy. 

Altogether, the total potential energy of the atomistic system is given by 



M-l 
j=-A/+l 

M 

i=-M+l 



M-l 

i=-Af+2 
2 



(3.5) 



yf - Oo 



ao 2 



We restrict ourselves to configurations in which the M leftmost atoms yf for — M-|-l<i<0 are 
situated in the interval {{i — f) 0,0, — ^) ao), whereas the M rightmost atoms yf for 1 < i < M 
are situated in the interval ((z — ^) ao, (« + ^) ao). The defect is situated between atoms y^ and 
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yf. In this case, the total energy simphfies to 

A/-1 M-1 

£^{yn = kf'i E ivhi - y? - ^0? + hk2 ivt+i - yt-i - ^^of 

i=-M+l i=-M+2 

M (3-6) 

i=-M+l i=l 

3.2. Quasicontinuum Approximation of the Frenkel-Kontorova Model. We now apply 
the quasicontinuum method to the dislocation model described in Section 13. 1[ The total energy 
(j3.6p is split up into atom-wise contributions, separately for the elastic interactions and the misfit 
interactions: 

^riy'') = ihiyt - yti - ao? + lk,{yt^i - vt " aof 

+ -My^ - VU - 2ao)' + \k2{y1^2 - yt - 2ao)', 

(3.7) 



an, we us( 
are neglected. We thus have that 



a.^c a-. ^ ) \kQ {yf -{i- l)aof , i = -M + 1, . . . , 0, 
' ' ^'ko{yf-iao)\ i = l,...,M. 



K 2 

To simplify notation, we use the convention that the undefined terms at the endpoints of the chain 



M 

f-(y«)=^-.™(y-)+£--.'=(y-)= ^ [^^'"(y'^) + (y"^)] . (3.8) 

i=-M+l 

Since the largest displacement of the atoms is to be expected near the defect, we deem the atoms 
—K + 1, . . . , K atomistic and the remaining atoms —M + 1, . . . , —K and K + 1, . . . , M continuum. 
Here K < M is some constant whose optimal value will be determined by the algorithm given in 
Section [5l 

The optimal choice of the repatoms for coarsening is investigated in the second paper of this 
series, so we work with a general formulation which holds for any values of ij for now. However, 
there are two restrictions on the coarsening. Since the atomistic region must not be coarsened and 
since we need full refinement in the vicinity of two atoms around the atomistic region due to the 
NNN interactions, we have that 

£j=j, j = -K-l,...,K + 2. (3.9) 

Second, we require that 

i-N+i = -M + 1, £_7v+2 = -M + 2, iN-i = M-l, and In = M (3.10) 

to incorporate the boundary conditions later. 

The elastic part S""'^ is translationally invariant, so we perform its QC approximation as described 
in the previous section. This leads to the continuum energy density 

ct>%r) = ^h{r - aof + \k2{2r - 2aof 

= ^ki2ir-aof (3.11) 

where ki2 := ki + 4/^2 • 

Regarding the misfit part £°''^, the above technique cannot be applied since the potential is 
not translationally invariant. However, there is a different summation technique to achieve a 
computationally efficient formulation which avoids the costly interpolation operator. 



ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 9 



To shorten the notation, we let indicate the sum in which the first term and the last term 
are only counted half: 



n—l 



^ ^ Xj . — 2'^1^ ~^ ^ ^ ■^i ~t~ 2"^** (3.12) 



where m < n and Xi € M. It is easy to verify that 

1:%^ = ^!^ and Y^^{m-^) = '^ (3.13) 

i=0 j=0 



for m > 0. 

For all pairs (j, j + 1) of continuum repatoms, we now reformulate all terms from (12. 7p which 
involve the interaction between and (-j+i- For j > 0, we get by definition (|2.4p of the operator 
L, by definition ([321) of and by ([31"^ that 

^J+l / p _ ■ ■ _ B \ 2 



I'^o Z^i'^j+i-V +2^0 ^;2 

, ,^ (yf - ^iQo)(y|+i - ^j+iao) ^'/,^ 

^ (;/f - ^jao)^ 2z.| + ^ (y|Ji-^+i«o)' 2z.| + 

2^0 -2 7 I" 2^0 r;2" 



i/^ 6 z^j' 6 

_ ^.„„\^,gc _/<.,,„„\„3 



A:, 



711 



6 
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For j < 0, we similarly obtain 



2 / > 



2/C0 



(y]^i-(£,+i-l)ao)2 2z.3 + i.^. 



(yf - (^,- - l)ao)(y|^i - (£,+1 - l)ao) z.^ - i^j 



6 



^(L-^'-^+i/y5^) = ^!'™(/y«^) and £'i^^{L^^^+^Iyi') = £^^'^,{Iy'i') , the QC approx- 
imation of the chain can be given by 

-K 



mce cjj^ { 



■ E 



' qc qc^ 

yj+1 - yj 



K 



N-1 
j=K 



(3.14) 



' qc 

yj+1 



+ €^{yT^y]li) 



Note that the interpolation ly^'^ does not have to be computed here since the relevant terms only 
depend on uncoarsened parts of the chain. 

Additionally, we will consider the atomistic-continuum approximation 

of the atomistic energy without coarsening. It is given exactly like the QC approximation (j3.14p 
with the only difference being that Vj = 1 and ij = j everywhere. 

4. Primal and Dual Problems 

4.1. Problem Setup. We are now ready to set up the problems we will solve. We are interested 
in finding the minimum of the energy ()3.14p subject to given boundary conditions. 

We give the boundary conditions by constraining the deformation of two atoms at each end of the 
chain. This guarantees that the potential with next-nearest-neighbor interactions can be directly 
applied to all non-boundary atoms without having to neglect interactions. We define the spaces 



ya 



l>2M 



d2M-4 



y«^:=M2^, Fo''^=R2^-^ (4.1) 

The spaces V"" and Vq will also be used for the uncoarsened atomistic-continuum potential S""^, 
so there is no need to define spaces V""*^ and Vq'^. We let y*'^ G V" denote any vector which has 
the desired boundary values y'^j^_^_i, y'^]^_^_2, Um-i^ ^^'^ Vm^ ™d we let y^'^'^ S V'^'^ by any vector 
satisfying (recall (jS.lOp ) 



bcq 



„bcq _ be 
y-N+2 ~ y-M+2' 



_ be 
Vn-i — Vm-i-: 

beq be 

vn ~ yM- 



(4.2) 
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For any vector y £ Vq, we denote the extension by zero boundary conditions to be Jy € V"', so 

Jy:=[0 y^ Oj^GM^A/^ ^43) 

and similarly we denote the extension by zero boundary conditions of y G Vq'^ to be J'^^y G V^"^. 
The spaces of admissible solutions are then given by JVq + y^'^ C V"' and J'^'^V^^ + y^'^'' C V'^'^, 
respectively. We note that : V"' — > Vq* is the restriction operator defined by 

iJ^y)j=yj fori = -M + 3,...,M-2. 

The minima y°, y'^'^, and y'^'^ of the energy functionals S"''^, and £'^^ given by (j3.6p . (|3.15p . 
and (13.14P subject to the above "clamped" boundary conditions are characterized as 

y'' := argmin £:"(y) G F", (4.4) 
y"^:= argmin ^"'=(y) G y^ (4.5) 
:= argmin S"^ (y) e . (4.6) 

We note that the minima are uniquely determined because E""^, and E'^'^ are strictly convex. 

4.2. Matrix Formulation. For the subsequent discussion, it will be convenient to reformulate the 
total energies in matrix notation: 

^'^(y) = ^(y - Si'')'^D'^^E<^D<^{y -s-) + \{y- h-^K^y - b'^), (4.7a) 

^"'^(y) = ^(y - s.^^fD'^^E'^'D^y - a'^) + i(y - h'^f K%y - b'^), (4.7b) 

f9C(y) = l(y _ ai'i^Y DlcT-^qcj^qc^^ _ ^qc^ ^ 1 _ ^Qcf j^qc^^ _ ^qcy ^4 7^) 

The matrices G m(2A^-i)x2M jjqc ^ ^{2N-i)x2N compute the distance between two adja- 
cent atomistic positions; the matrices E" G m(2A/-i)x(2M-i)^ ^ac g ^{2M-i)x{2M-i) ^ ^qc ^ 
^{2N-i)x{2N-i) contain the spring constants h, k2, and k^; and the matrices iT" G K(2M-i)x(2Af-i) 
and G m(2^-i)x(2^-i) contain the misfit constant ko. The vectors a",b" G M^^^ and a«^ b^'^ G 
are constants describing the minimum energy deformations for the elastic energy and mis- 
fit energy. The precise and lengthy definitions for all of these matrices and vectors are given in 
Appendix [Al 

If we decompose y" = Jy"' + y^^ for y" G Vq, then the minimization problem is given as 

y'^ = argmin^"(Jy + y''^) 
yeVff 



arg mm 

yev-o" 



\{Jy + y^^- a") D^^E^D" {Jy + y^'^ - a'^) 

+ Wy + y'" - ^"fK"{jy + y^^ - b")l . (4.8) 



We also decompose y"^ = Jy"" + y^'^ and y"^ = Ji^y^"^ + y^^'' for y"" G and y''^ G V^^, and 
we then formulate similar minimization problems for y"'^ and y'^'^. Therefore, y*^, y"'^, and y'^" are 
determined by the linear systems 

M"y" = r, (4.9a) 
Miy}" = fi", (4.9c) 
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where 



J^iD^^E^D" + K'')J, 

-J^D'^^E'^^D'^iy'"' - a") - J^K°(y'"= - b"), 

_ JQCT jjqcT j^qc jjqc ^ybcq _ ^qc-^ _ jqcT j^qc ^^bcq _ y^qc^ 



(4.10) 



We note that the matrices M", M"''^, and M'^'^ are positive definite, so the total energies admit a 
single global minimum and no other local minimum. 

4.3. Goal-Oriented Error Estimation. To compare the approximate QC model to the original 
atomistic model, we have to analyze how much the solution of the atomistic model deviates 
from the solution y"?*^ of the QC model. This deviation, which can be viewed as an approximation 
error, can be measured in different ways, for example as ||y" — IJ'^'^y'^'^W for some norm ||-||. Here 
we follow a different approach, namely we measure the error of a quantity of interest denoted by 
Q(y) for some function Q : ]R^*^~^ — > M. Hence, we intend to estimate 

Q(y") - Q(J^/J''V)- (4.11) 

We will assume for simplicity that Q is linear and thus has a representation Q{y) = q^y for some 
vector q G Vq. 

For our application, a natural quantity of interest is the size of the dislocation, that is, the 
distance between the two atoms yo and yi to the left and right of the dislocation. This gives us 

Q(y) = q^y = yi - yo with q = [o, . . . ,o, -i, 1,0, . . . ,o]^. (4.12) 

Two different sources of error arise during the QC approximation, namely the localization of the 
potential energy, that is, the passage from the atomistic to the continuum formulation on the one 
hand, and the coarsening in the continuum region by the restriction to the repatoms on the other 
hand. We denote these two errors by 

e := y" - y'^^ and e""'^'' := y"" - J^/J^'^y'?^ (4.13) 

It makes sense to study these sources independently. Employing the linearity of Q, we have that 

\Q{yn - QiJ^Wym = IQ(e) + Qie'^'^ni < IQ(e)| + IQie'^'ni- (4.14) 
The error term |(5(6)l will be studied in Section [5l and the error term |Q(e°^*'^)| will be studied in 
the second part of this paper series. 

4.4. Dual Problems. To facilitate the goal-oriented error analysis, we introduce the dual problems 



M^g"" = q, (4.15a) 
^jacgac ^ ^ (4.15b) 

for g'^,g"'^ E M^^'^~^, and g"^"^ G M^^"^. We note that the dual problems differ from the primal 
problems only by the right hand side since the matrices M", M""^, and M'^^ are symmetric. 

The solutions g", g"'^ and g'^^ can be viewed as influence functions: They describe how the error 
at a specific point in the domain influences the error measured in terms of the goal function. 
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Analogously to the primal errors (|4.13p . we define the dual errors 

e := g« - g'^^ and e"^"?^ := g"^ - J^IJ^^g'^". (4.16) 
In addition, we will need the primal and dual residuals 



i?"(g) 
i?-(g) 



(y"^ - y) = - My, 
M^" (y"^ - y) = - iW^^y, 
(g« - g) = q _ M«g, 



5. Error Estimation for Atomistic vs. Continuum Modeling 

In this section, we estimate the error arising from the approximation of an atomistic model 

by a continuum model. We consider y'*'^ and g"'^ to be computable, although in practice we can 
only compute the coarsened approximations y"^"^ and g'^^. 

To this end, we adapt a technique introduced in [l^ and [l^ to estimate the modeling error for 
an elasticity model with rapidly oscillating coefficients and its homogenized version. We generalize 
this technique such that it allows for different right hand sides f° and f"'^ of the primal problem 
()4.9p instead of a common right hand side as it is used in the above-mentioned works. 

We have 

^(y") - Qiy"! = q^e = g^'^M^'e = (g^^^ + e)M"e 

= g"^^i?"(y''^) + e^We. (5.1) 

The term g^cT^aj^yac^ 

can be computed, whereas e-^M"e cannot because both e and e are nu- 
merically unknown. Instead, we estimate e^M"'e from above and from below by quantities that 
actually can be computed. 

We will give two different error estimators rji and r]2- Before, we need to derive some auxiliary 
estimates to facilitate their development and analysis. 

5.1. Auxiliary Estimates. We reformulate the difference y"^ — y'^'^ of the respective solutions in 
terms of a difference of the energy matrices. To this end, we define the perturbation matrix 

P :=! - {E'')-^E'"' (5.2) 

where I denotes the identity matrix. Note that EaP = E"" — E""^. 

Lemma 5.1. For any a, /3 G M, we have that 

M'^(ae + /3e) = -J^D^^E^PD" [a( Jy'^^ + y^'^ - a'^) + ^Jg""] . (5.3) 

Proof. We conclude from (j4.9p that 

M^e = My - M'^'y'"' + {M'"' - M") y""" 

and similarly since M"g" = M^'^g"'^ = q that 

M"e = M°(g'^ - g"^) = (iVr^ - iVr)g'^^ (5.5) 
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Thus, it follows from (flTO|) and ([O]) that 

M^'iae + pe) = a [{M'"' - M'')y''^ + P - + PiM"" - M'')g'"' 

= J^D^^iE"^ - E'')D'' [a( Jy''" + y'"^ - a'') + /3Jg"'=] (5.6) 

We note that the iC'^-related terms cancel here, because they coincide for the atomistic model and 
the continuum model. □ 

Lemma 5.2. We have that 

\\ae + Pe\\M^ < ||PZ)'^[a(Jy'^^ + y"^ - a'^) + /JJg'^^] ||^,. (5.7) 

We note that the right hand side is numerically computable. 

Proof. To shorten the notation, we abbreviate z = a(Jy'^^ + y*'^ — a") + pjg"''^. By Lemma lS-H we 
have 

v^M'^(ae + /?e) 
||ae + pe||M<^ = sup 1| — 1| 

= sup — 

veVo^Vlo} IKI|m« (5.8) 

\\D''Jv\\Ea\\PD''z\\Ea 

— iTf; — r~n 

veV(f\{0} ||iJ"Jv||£;a 

= WPD^'zWEa. 

Here we have used that HZ?" Jv||£;a < ||v||Afa because the matrix K"^ in ()4.10p is positive definite. □ 

5.2. First Error Estimator. We are now ready to derive the first error estimator, iji. By the 
parallelogram identity, we have for all cr that 

e^M"e = (CT-^e^)M"((je) 

II / -1-112 HI -1NI2 (5.9) 

In the following, we will determine computable constants r]^^, t]^^, rj^pp and 7]~pp such that 

??iow ^ Ike + cr"^e||A/« < r?+pp' 
< Ike - o-~^e||A/" < Vnpp- 
From Lemma 15.21 we immediately get the upper estimates rj^pp and 77^pp : 



(5.10) 



??+ p := \\PD'' [a( Jy-^^ + y^'^ - a'^) + a" Vg'^^] 
??-pp := ||PD" [a{Jy-' + y'"^ - a'^) - a'^Jg^^' 



(5.11) 



We note that ^io„, ??upp and r/~pp will depend on cr, but the estimates (|5.10p will hold for any 
o" 7^ 0. We will now choose a in such a way that the estimates are as sharp as possible, that is, 
such that ry+pp and r/~pp are smallest. 

Lemma 5.3. Both ij^pp and ij^^p given by (j5.1ip attain their minima for 



_ / llPD^'Jg^'^lls" 
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Proof. We have that 

(7y± p)2 = a^\\PD''{Jy'''' + y^'^ - a'^)|||„ ± 2g'"''^ j'^ D^'^ E"" P D" {Jy""^ + y''^ - a'^) 

+ a-2||PZ)Vg'^i|,. (5.13) 

Setting the first derivative of the mapping a i— > (?/upp)^ to zero, we obtain the condition 

2a||PZ5'^(Jy"^ + y''^-a")|||, -2a-3||PZ)Vg"i|„ = (5.14) 

for critical points of (??upp)^- This equation has the unique positive solution ()5.12p . Because 
lim|o-|_^oo ry^pp = linio-^o ??^pp = oo, this point corresponds to a minimum. Hence the quantities ry^pp 
attain their minima at cr = ct. □ 



Regarding the lower bounds rj^^ and ^i^^, we have 

„ v^M"(CTe ± CT-^e) 

llcreiba" e\\M^ = sup — 7 

vev;f\{o} 

v^(ai?"(y-) ± a-ii?«(g-)) ^ 

= sup 7 (5.15) 

^ v^(CTi?'^(y"'=) ± a-ii?"(g'^^)) 
~ llvolUf'^ 

for any vector vq G Vq \ {0}. Numerically, we have the two vectors y'^'^ and g""^ at our disposal, 
hence it makes sense to take a linear combination vq = y"'^ + O^g,"''^. Here we follow the strategy 
of 

and choose 6^ as the critical points of Vi^^- 



Lemma 5.4. Let 



Then the lower bounds 



have a unique critical point for 



aR^iy"^) ± a-^R^g""). (5.16) 



_itT-.ac „acT ]\/fa^,ac „±r„ac ||--ac||2 

j± r y g y — r g ||y 



(5.18) 



Proof. We have 

d^''^"- " !|y- + ^±g-||i,„ ■ ^^"^^^ 

Setting this expression to zero and solving for 9^ leads to the above condition. □ 

However, let us note that this critical point is not necessarily a maximum of rj^^, which would 
be optimal for bound (|5.10p . Depending on the actual vectors y""^ and g"'^, it can be shown that 
this critical point could be a minimum. 

Now we have all necessary ingredients to construct the error estimator rji. From (|5.ip and (j5.9p . 
we get the computable estimate 

^acTj^a^yac^ + U^L? " i(^upp)' < ^?(y") " Q(y"') < g"^^i?'^(y"^) + l{vtr>p f " K^l^w)'' 

(5.20) 
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At first sight, this looks hke we could get an estimate for |Q(y") — Q(y"'^)| from both above and 
below. However, this is only true if both the left hand side and the right hand side have the same 
sign, which in general does not hold. But we get the following estimate. 

Theorem 5.1. We have that 

\Q{yl-Q{yn\<vi, (5.2i) 

where the computable error estimator is defined as 

m := max (|g-T^a(yac) + 1(^+^)2 _ 1 (^-^^)2| ^ IgacT^a(yac) + 1(^+^)2 _ 1 (^-j2|) . (5.22) 

We note that the computation of the rj^pp terms involves the solution of a linear system with 
matrix E"" as its inverse appears in the operator P. The matrix is not diagonal, but has 
condition number So this is negligible compared to what would be necessary to solve the 

original atomistic problem which includes the operator D""^ E°'D°' with condition number 0{M'^). 



5.3. Second Error Estimator. There is no reasonable way to decompose the error estimator 771 
into a sum of element-wise or atom-wise contributions due to the f]^^ terms. Therefore, we derive 
another error estimator 772 which allows for such a decomposition, at the price of a less accurate 
estimate than 771. 

Theorem 5.2. We have that 



M-2 



\Q{y')-Q{yn\<V2< 

i=-A/+3 



V2A 



M-1 

+ E 

i=-M+l 



V2,i 



(5.23) 



where the computable global error estimator r]2 and the computable local error estimators, and 
r]f-, associated with atoms and elements, respectively, are defined as 



m 



jgacT^a^yac^l _^ \\P {Jy"" + y''" - ai'')\\Ea\\P Jg'^'^WEa , 
|^ac^a(yac).| ^ Z = -M + 3, . . . , M - 2, 

i I (PL>"(Jy"^ + y^'^ - a")) .((^" - £;"^)D''( Jy"^ + y''^ - a'')) .| 
+i |(PZ)"Jg'^^)i((^'^ - E'"')D''Jg'"').\ , i = -M + 1,...,M -1. 



Proof. From (j5.ip and Lemma 15.21 we conclude that 



(5.24) 



IQiy") - Qiy"')] < |g'''^^''(y"')| + |e^Are| 

< |g"'=^i?"(y"^)| + ||e||Af"||e||M« 

< |g'''^i?"(y''')| + WPD^iJy"" + y^'' 



lJa\\PD''Jg'"'\\E'^ 



(5.25) 
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(5.26) 



which gives us the global estimate. For the decomposition into local contributions, we further 
estimate 

< |g'^'=^i?"(y°")| + i||PZ)'^(Jy"^ + y^= - a'^)|l|„ + ijl Pi? Vg'^^|||. 

M-2 M-1 

< Yl larR'iynil + h E |(^'^Vg'^^)i((i^"-E'^^)DVg"^).| 

i=-A/+3 i=-M+l 
M-l 

+ i |(^^"('/y"^ + y*''-a")).((^''-^''')i:'"(Jy"' + y^'=-a"))^ 

j=-A/+l 
M-2 M-l 

i=-M+3 i=-M+l 

which completes the proof. □ 

Let us remark that instead of the first inequality in (j5.26p . one can get an apparently better 
estimate 

\\PD''{Jy''^ + y''^ - a'^)||sa WPD" J^'^'Wec^ 

= \^\\PD-{Jy-- + y"^ - a'^)!!!. + \^-^\\PD- J^^'Wl^ ^^'^^^ 
by introducing the additional weight factor 

llPP'^Jg'^^lliJa 

and then decomposing the resulting terms similar to the above. However, our numerical results 
showed that this modification does not significantly improve the decomposed error estimator for 
the application considered here. 

6. Numerics 

In the preceding sections, we constructed the error estimators rji and r]2. We will now give an 
algorithm for adaptive atomistic-continuum modeling based on these error estimators. Then we 
will present and discuss some numerical results. 

6.1. Algorithm. The error estimator r]i should give a better estimate of the error than t/2, because 
r/2 involves the inequality |e^M"e| < ||e||Afa ||e||M" in (|5.25p in contrast to the parallelogram identity 
for r]i. However, r]2 can be decomposed into atom- wise and element- wise contributions r^g* and 
rj2i, whereas the r]^^ terms in rji do not admit a reasonable decomposition that can be used for 
atomistic-continuum adaptivity. 

We make use of this by employing the sharper estimate rji to determine whether a given global 
error tolerance Tgi for the error in an adaptive algorithm has already been achieved or not. If 
not, we use the decomposed estimates r/f* and rjfi to determine where the more precise atomistic 
modeling is needed. This leads us to the following algorithm: 

(1) Choose Tgi. Model all atoms as a continuum. Set Tat ^ Tgi- 

(2) Solve primal problem ()4.9bp for y"'^ and dual problem ()4.15bp for g^"^. 

(3) Compute error estimator rji from (j5.22p . 

(4) If rji < Tgi, then stop. 
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M 


iteration 


K 


Tat 


'^^ 

1 J- 


100 


1 





l.OOOOOOe-10 


3.899207e-02 




2 


28 


l.OOOOOOe-11 


5.915080e-10 




3 


32 


l.OOOOOOe-12 


4.878532e-ll 


1000 


1 





l.OOOOOOe-10 


3.899208e-02 




2 


28 


l.OOOOOOe-11 


5.915100e-10 




3 


32 


l.OOOOOOe-12 


4.878548e-ll 


10000 


1 





l.OOOOOOe-10 


3.899208e-02 




2 


28 


l.OOOOOOe-11 


5.915100e-10 




3 


32 


l.OOOOOOe-12 


4.878548e-ll 


100000 


1 





l.OOOOOOe-10 


3.899208e-02 




2 


28 


l.OOOOOOe-11 


5.915099e-10 




3 


32 


l.OOOOOOe-12 


4.878540e-ll 


1000000 


1 





l.OOOOOOe-10 


3.899208e-02 




2 


28 


l.OOOOOOe-11 


5.914422e-10 




3 


32 


l.OOOOOOe-12 


4.871775e-ll 



Table 1: Convergence of the algorithm for Tgi — 10 and different values of M. 



(5) Compute local error estimators "q^,- and r/|'^ from (j5.24p . 

(6) Set Tat ^ 

(7) Make all atoms i atomistic for which 

Vf^ ■■= ri% + \ {vlU + Vti) > Tat- (6.1) 

(8) Go to (2). 

Here T^iv > 1 is a constant factor which describes how fast the atom-wise tolerance Tat should 
decrease during adaption. Our experience has been that t^h^j = 10 is a reasonable choice. 

The crucial adaption step is (7). The adaption criterion ()6.ip deems all atoms to be modeled 
atomistically if the associated error from the decomposition of ri2 exceeds the atomistic error toler- 
ance Tat- Here the element-wise errors rjf^ are distributed equally to the two adjacent atoms i and 
i + l. 

For the dislocation at the center of the chain and the chosen goal function, we expect that the 
atomistic repatoms always form a symmetric interval around the center. We have used the above 
adaptive atomistic-continuum algorithm to approximate our Frenkel-Kontorova model and have 
always found that the atomistic region is the set of atoms —K + 1,...,K for some K depending 
on M and Tgi. Thus, the modeling approach given in Section [3] of restricting to an atomistic region 
consisting of atoms —K + 1, . . . ,K for some K rather than considering a more general atomistic 
region is justified a posteriori. 

6.2. Numerical Results. The algorithm has been implemented as described above. The boundary 
conditions were chosen as 

y'L'M+i = -M, y%,+2 = -M + l, 2/^^_i = M-l, y^ = M. (6.2) 

The elastic constants are ko = 1 and ki = k2 = 2. 



10"' 



10 



10" 



10"' 



10 ■ 
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10" 




10 



10" 



10" 



10 




-500 -400 -300 -200 -100 100 200 300 400 500 -500 -400 -300 -200 -100 100 200 300 400 500 

Figure 4: The error estimators 773°* (left) and r?!',, r;^*, (right) for M = 500, K ^ 20. 



Table [T] shows how the algorithm given above performs. After 3 iterations, the desired accuracy 
Tgi = 10"^'^ is achieved. Moreover, we can see from the table that the number of iterations are 
independent of M, that means the algorithm behaves robustly with respect to the problem size M. 

Figure H] (left) shows the decomposition of the error estimator rj2 for a typical setting M = 500, 
K = 20. One can clearly see that the error in the atomistic region is small, whereas the error is 
large in the continuum regions that border the atomistic region. It then decreases exponentially 
towards the endpoints. 

The error in both the atomistic region around the center and the continuum regions far away 
from the center are in the range of the (relative) machine precision Emach, which accounts for the 
fluctuations in these regions. The error can be considered to be numerically zero in these regions. 
In the continuum regions, we observe an error of magnitude ^(Smac/i)' whereas in the continuum 
region we have O (emach), which leads to the different magnitudes of the fluctuations. 

Figure S] (right) shows the element-wise contributions r/^'^ and the atom-wise contributions r/2* 
of the decomposed error estimator 7/2°* = r/f* + ^{vfi-i + vfi)- "^^^ atomistic part r/2*, which 
corresponds to the ^"■'^'^ R°- {y"-'^) term, is dominant in the sense that it is about ten times larger 
than r/2'^, which comes from the estimate for the perturbation term e^M"e. The fluctuations 
due to the limited machine precion in the atomistic region come from rj^^^ whereas those in the 
continuum region away from the defect stem from 77!'^. Let us note that in other applications of 
duality-based error estimation, the first term might not always be the dominant term. For example, 
in mesh refinement for classical linear finite elements, the first term even vanishes due to Galerkin 
orthogonality. 

Table [2] and Figure [5l which display the same data, show the efficiency of the error estimators, 
^i/IQ(y" ~ y"'^)! ^-iid 7?2/|Q(y" — y"'^)!) for M = lOOO. For comparison, the actual error is given as 
well. For the relatively small ID problem, the actual error can be easily computed, whereas in real 
world applications it is of course not available. One can clearly see that 771 gives a better estimate 
than ?72, which numerically confirms our conjecture that rji is a better estimator than rj2- We see 
that r/i overestimates the actual error by a factor of 1.4, while 772 is in a still acceptable range of 
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K 




1C\ 

/ 




'11/ J / 


V2 









3.627633e- 


02 


3 


899208e- 


02 


1.074863 


3 


999783e- 


02 


1 102588 


2 


3.375762e- 


02 


3 


872272e- 


02 


1.147081 


5 


1017006- 


02 


1.511274 


4 


3.468605e- 


03 


4 


3435956- 


03 


1.252260 


5 


422007e- 


03 


1 563166 


6 


5.418585e- 


04 


7 


156249e- 


04 


1.320686 


9 


187940e- 


04 


1 695635 


8 


1.227067e- 


04 


1 


675383e- 


04 


1.365356 


2 


193196e- 


04 


1.787348 


10 


3.287188e- 


05 


4 


540984e- 


05 


1.381419 


5 


984186e- 


05 


1.820457 


15 


1.416914e- 


06 


1 


966114e- 


06 


1.387603 


2 


597488e- 


06 


1.833201 


20 


6.267636e- 


08 


8 


695824e- 


08 


1.387417 


1 


148736e- 


07 


1.832805 


25 


2.770161e- 


09 


3 


843388e- 


09 


1.387424 


5 


077204e- 


09 


1.832819 


30 


1.224369e- 


10 


1 


698739e- 


10 


1.387440 


2 


244073e- 


10 


1.832840 


35 


5.410783e- 


12 


7 


508365e- 


12 


1.387667 


9 


918687e- 


12 


1.833133 


40 


2.379208e- 


13 


3 


318024e- 


13 


1.394592 


4 


383361e- 


13 


1.842362 


45 


8.992806e- 


15 


1 


430733e- 


14 


1.590975 


1 


901601e- 


14 


2.114580 


50 


7.771561e- 


16 


4 


120094e- 


16 


0.530150 


6 


201285e- 


16 


0.797946 



Table 2: Efficiency of the error estimators, r?i/|(3(y'' - y'''^)! and r72/|Q(y° - y""")!, for M = 1000. For 
if = 45 and K = the results become inaccurate due to hmited machine precision. 
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Figure 5: Efficiency of the error estimators for M = 1000. 
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optimal K 


K by ?7i 


K bv 779 


le-02 


3 


3 


3 


le-03 


5 


5 


5 


le-04 


9 


9 


10 


le-05 


12 


13 


13 


le-06 


16 


17 


17 


le-07 


20 


20 


21 


le-08 


23 


24 


24 


le-09 


27 


28 


28 


le-10 


31 


31 


32 


le-11 


35 


35 


35 


le-12 


38 


39 


39 


le-13 


42 


42 


43 


le-14 


45 


46 


47 



Table 3: Efficiency of the error estimators for M — 1000. 



2 times the actual error. Moreover, we can see from Table [2] and Figure [5] that the error decreases 
exponentially with K. 

Finally, we compare the optimal (smallest) value of K which is needed to achieve a given accuracy 
Tgi with the values for K determined by the error estimators rji and rj2, again taking into account 
the precise error which is available for the model problem. We see from Table El that even t]2 only 
overestimates K by at most 2 atoms. Thus, we get an efficient estimate of the required atomistic 
region for both error estimators. 

Appendix A. Matrix Definitions 



We describe the matrices from Section 14.21 The matrix 

-1 1 



-1 1 



-1 1 



o{2M-l)x2M 



(A.l) 



transforms atomistic positions to distances between adjacent atoms. Similarly, 

^ ^{2N-l)x2N 



,-1 
-N+1 



,-1 

"^-Ar+2 '^-N+2 



-1 



-1 



(A.2) 



transforms repatom positions from a coarsened chain to normalized distances between adjacent 
repatoms. The matrices 

'ki + k2 i= j e {-M + l,M -1} 
ki + 2k2 i = j e {-M + 2,...,M -2} 
h \j -i\ = l 

otherwise. 



(A.3) 
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and 



qc\ 







ifc2(^?_l + <5f+l) 





= < 



j = i + l 
j = i-l 
otherwise, 

j 

-i-l 



(A.4) 



I ■ 

j 
j 

otherwise, 



(A.5) 



for i,j = —M + 1, . . . , M — 1 and i,j = —N + 1, . . . ,N — 1, respectively, describe the spring 
interactions in terms of the distances between atoms or repatoms. Accordingly, the matrices 

ko 



ko 



and 



^Iko [i2l^^.l + u-_\) + {2l,, + l^r^)] 
lko{2u_N+i + i^Zn+i) 
lko{2uN-i + z^j^ii) 
^ko{i^i - 
^ko{i^j - 1^7^) 




r,2Mx2M 



i = jG{-N + 2,...,N-l} 
i=j= -N + 1 
i = j = N 
j = i + l 
j = 1-1 



(A.6) 



(A.7) 



otherwise, 



for i, J = —N + 1, . . . , describe the misfit interactions for the original atomistic system and the 
QC approximation. Finally, the constant vectors 



a'* = [(-M + l)ao (-M + 2)oo • • • (M - l)ao Mao] ^ e 



l)2M 



D2Ar 



d2M 



(A.8a) 
(A.8b) 
(A.8c) 

^2M 

(A.'sd) 



b" = [-Mao (-M + l)ao • • • -oq Oq • • • (M - l)ao Moq]^ G 
b^^ = [(£_Ar+i - l)ao (■^-Ar+2 - l)ao ••• (4 - l)ao ^1^0 •■■ -^Ar-iao ^Araoj^el 

fix the equilibrium positions for the spring interactions and the misfit energy, respectively. 
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